function [VARcoeff , VARcoeff_adj , VARerror_adj] = VAR1_bsadj(data , B , method , trend)
% Estiamte a VAR(1) for input data and use bootstrap to bias-adjust the
% coefficients.
% Require James P. LeSage's Econometrics Toolbox
% INPUTS:
%   DATA   - T by K vector of data to be bootstrap adjusted  
%   B      - Number of simulations
%   method - Method of bootstrap data residuals
%          - 1: i.i.d. bootstrap
%          - 2: circular bootstrap
%          - 3: stationary bootstrap
%   trend  - whether include a trend in estimating VAR[OPTIONAL]. Default
%            no trend: 0. To include a trend, put trend = 1.
% 
% OUTPUTS:
%   VARcoeff     - OLS estimates of VAR coefficients
%   VARcoeff_adj - Bootstrapped bias-adjusted VAR coefficients
%   VARerror_adj - Bootstrapped bias-adjusted VAR residuals
% 
%% ### Input Checking

if nargin < 3
    error('3 inputs required')
end

% Get the length of the data
[T,K]=size(data);

if T < 2 
    error('DATA must have at least 2 observations.')
end

switch nargin
    case 3
        trend = 0;
end

%% ### For later use

% RHS constant and trend term
if trend == 1
    RHS_exo = [ones(T - 1 , 1) , (2 : T)'];
elseif trend == 0 
    RHS_exo = ones(T - 1 , 1);
end

% Optimal block length
L_optimal = opt_block_length(data); %See Andrew Patton's Website
  
% Block length. Following Crump and Gospodinov (2019) to use max 
if method == 1 % i.i.d
   L = 1;
   type = 1;
elseif method == 2 % circular
   L = ceil(max(L_optimal(2 , :)));
   type = 2;
elseif method == 3 % stationary
   L = ceil(max(L_optimal(1 , :)));
   type = 3;
end

%% Estimate VAR 
% Raw estimates using OLS
b_raw = olsgmm(data(2 : end , :) , [RHS_exo , data(1 : end - 1 , :)] , 0 , -1);

e_raw = data(2 : end , :)...
      - [RHS_exo , data(1 : end - 1 , :)] * b_raw;
  
  
%% Finite sample adjustment

b_bs = nan(K * (1 + trend) + K^2 , B);

% Bootstrap errors (index)
error_index = bs_choose((1 : T - 1)' , B , L , type);

% Create bootstrap samples
% Following Kothari and Shanken(1997), condition on the first observation
% We do not directly reorder element but reorder errors to preserve trend
for i = 1 : B    
    e_bs_temp = e_raw(error_index(: , i) , :);
    
    data_bs_temp = nan(T , K);
    data_bs_temp(1 , :) = data(1 , :);
    for t = 2 : T
        data_bs_temp(t , :) = [RHS_exo(t - 1 , :) , data_bs_temp(t - 1 , :)] * b_raw ...
                            + e_bs_temp(t - 1 , :);
    end    
    % Re-estimate VAR coefficients
    b_bs_temp = olsgmm(data_bs_temp(2 : end , :) , [RHS_exo , data_bs_temp(1 : end - 1 , :)] , 0 , -1);
    
    % store coefficients
    b_bs(: , i) = b_bs_temp(:);
end

% E(\hat{t}) - t_0 = E(\tilde{t}) - \hat{t}
% t_0 in vector

b_adj_vec = 2 * b_raw(:) - mean(b_bs , 2); % Bias

VARcoeff = b_raw;

VARcoeff_adj = reshape(b_adj_vec , [] , K);

VARerror_adj = data(2 : end , :)...
              - [RHS_exo , data(1 : end - 1 , :)] * VARcoeff_adj;
          
end
